global Y T1 T2 SM SR Draws M Xc Gmax Gmin ll or orm

% 
% wparam = BetaProbitWhite;
bparam = BetaProbitBlack(1:27);
% 
% Sw = SigmaProbitWhite;
% 
% theta2 = [BetaProbitBlack;0;0;0;0;0;0];
% I=(cbyrace3==1)+0;
% 
% Y=Yo(I==1);
% T1=T1o(I==1);
% T2=T2o(I==1);
% SM=SMo(I==1);
% SR=SRo(I==1);
% Draws = Drawso(I==1,:);
% Xc=Xco(I==1);
% or = repmat((cbyerace3(I==1)==0)+0,1,M);
% orm = repmat((cbymrace3(I==1)==0)+0,1,M);
% 
% 
% Sb=ERRORS(BetaProbitBlack(1:27));
% 
% 
% theta = Drawso;
% Ib = (cbyrace3==1)+0;
% Iw = (cbyrace6==1)+0;
% T1=T1o;T2=T2o;G=Xco;
% 
% bias1 = repmat(T1,1,500)-normcdf((bparam(20)+bparam(21)*repmat(G,1,500)).*repmat(Ib,1,500)+(wparam(20)+wparam(21)*repmat(G,1,500)).*repmat(Iw,1,500));
% bias2 = repmat(T2,1,500)-normcdf((bparam(20)+bparam(21)*repmat(G,1,500)).*repmat(Ib,1,500)+(wparam(20)+wparam(21)*repmat(G,1,500)).*repmat(Iw,1,500));



% cm = param(1);
% cr = param(2);
% phim = param(3);
% phir = param(4);
% sigem = param(5);
% siger = param(6);
% cgpa = param(7);
% phigpa = param(8);
% sigegpa = param(9);
% 
% c1 = param(10);
% c2 = param(11);
% b1 = param(12);
% b2 = param(13);
% phi1 = param(14);
% phi2 = param(15);
% gamma1 = param(16);
% gamma2 = param(17);
% sigtheta = param(18);
% sige12 = param(19);
% cy = param(20);
% by = param(21);
% 
% sige1 = 1;
% sige2 = 1;
% sigy = 1;
% 
% c1d = param(22);
% c2d = param(23);
% phi1d = param(24);
% phi2d = param(25);
% b1d = param(26);
% b2d = param(27);

% 
% for i = 1:2
%     if i==1;
%         I=cbyrace6==1;
%         G=Xco(I==1,:);T1=T1o(I==1,:);T2=T2o(I==1,:); SIG = Sw;b1=bias1(I==1,:); b2=bias2(I==1,:);
%         param=wparam;b1=mean(mean(bias1)); b2=mean(mean(bias2));
%     elseif i==2;
%         I=cbyrace3==1;
%         G=Xco(I==1,:);T1=T1o(I==1,:);T2=T2o(I==1,:); SIG = Sb;b1=bias1(I==1,:); b2=bias2(I==1,:);
%         param=bparam;b1=mean(mean(bias1)); b2=mean(mean(bias2));
%     end
%     
% %     dYdb1(i) = mean(mean(param(13)*normpdf(param(11)+param(13)*b1+param(14)*b2+param(15)*repmat(G,1,500))));
% %     dYdb2(i) = mean(mean(param(14)*normpdf(param(11)+param(13)*b1+param(14)*b2+param(15)*repmat(G,1,500))));
% 
%     
%     
%     dYdb1(i) = mean(param(16)*normpdf(param(20)+param(16)*b1+param(17)*b2+param(21)*G));
%     dYdb2(i) = mean(param(17)*normpdf(param(20)+param(16)*b1+param(17)*b2+param(21)*G));
%     Jyt1 = NumJacobian(@(par)...
%         mean(mean(par(16)*normpdf(par(20)+par(16)*b1+par(17)*b2+par(21)*G))),param);
%     Jyt2 = NumJacobian(@(par)...
%         mean(mean(par(17)*normpdf(par(20)+par(16)*b1+par(17)*b2+par(21)*G))),param);
% 
%     
%     
% %     Jyt1 = NumJacobian(@(par)...
% %         mean(mean(par(13)*normpdf(par(11)+par(13)*b1+par(14)*b2+par(15)*repmat(G,1,500)))),param);
% %     Jyt2 = NumJacobian(@(par)...
% %         mean(mean(par(14)*normpdf(par(11)+par(13)*b1+par(14)*b2+par(15)*repmat(G,1,500)))),param);
%     Vyt1o(i) = sqrt(Jyt1*SIG*Jyt1');
%     Vyt2o(i) = sqrt(Jyt2*SIG*Jyt2');        
% end
% 
% [dYdb1;Vyt1o;dYdb2;Vyt2o]




for i=1:2;
    if i==1;
        I=cbyrace6==1;
        G=Xco(I==1,:);T1=T1o(I==1,:);T2=T2o(I==1,:); SIG = Sw;b1=bias1(I==1,:); b2=bias2(I==1,:);
        param=wparam;b1=mean(mean(bias1)); b2=mean(mean(bias2));
    elseif i==2;
        I=cbyrace3==1;
        G=Xco(I==1,:);T1=T1o(I==1,:);T2=T2o(I==1,:); SIG = Sb;b1=bias1(I==1,:); b2=bias2(I==1,:);
        param=bparam;b1=mean(mean(bias1)); b2=mean(mean(bias2));
    end
    ela1(i) = mean((param(16)*b1)*normpdf(param(20)+G*param(21)+param(16)*b1+param(17)*b2)./normcdf(param(20)+G*param(21)+param(16)*b1+param(17)*b2));
    ela2(i) = mean((param(17)*b2)*normpdf(param(20)+G*param(21)+param(16)*b1+param(17)*b2)./normcdf(param(20)+G*param(21)+param(16)*b1+param(17)*b2));
    
    
    Jyt1 = NumJacobian(@(par)...
        mean((par(16)*b1)*normpdf(par(20)+G*par(21)+par(16)*b1+par(17)*b2)./normcdf(par(20)+G*par(21)+par(16)*b1+par(17)*b2)),param);

    Jyt2 = NumJacobian(@(par)...
        mean((par(17)*b2)*normpdf(par(20)+G*par(21)+par(16)*b1+par(17)*b2)./normcdf(par(20)+G*par(21)+par(16)*b1+par(17)*b2)),param);
    
    

    
    
    Vyt1o(i) = sqrt(Jyt1*SIG*Jyt1');
    Vyt2o(i) = sqrt(Jyt2*SIG*Jyt2');        
    
        
end

[ela1;Vyt1o;ela2;Vyt2o]


% cm = param(1);
% cr = param(2);
% phim = param(3);
% phir = param(4);
% sigem = param(5);
% siger = param(6);
% cgpa = param(7);
% phigpa = param(8);
% sigegpa = param(9);
% 
% c1 = param(10);
% c2 = param(11);
% b1 = param(12);
% b2 = param(13);
% phi1 = param(14);
% phi2 = param(15);
% gamma1 = param(16);
% gamma2 = param(17);
% sigtheta = param(18);
% sige12 = param(19);
% cy = param(20);
% by = param(21);
% 
% sige1 = 1;
% sige2 = 1;
% sigy = 1;
% 
% c1d = param(22);
% c2d = param(23);
% phi1d = param(24);
% phi2d = param(25);
% b1d = param(26);
% b2d = param(27);
    

theta1=[wparam(20);bparam(20)];
theta1=[0;0];

for i=1:2;
    if i==1;
        I=cbyrace6==1;
        G=Xco(I==1,:);T1=T1o(I==1,:);T2=T2o(I==1,:); SIG = Sw;b1=bias1(I==1,:); b2=bias2(I==1,:);
        param=wparam;b1=mean(mean(bias1)); b2=mean(mean(bias2));
    elseif i==2;
        I=cbyrace3==1;
        G=Xco(I==1,:);T1=T1o(I==1,:);T2=T2o(I==1,:); SIG = Sb;b1=bias1(I==1,:); b2=bias2(I==1,:);
        param=bparam;b1=mean(mean(bias1)); b2=mean(mean(bias2));
    end
    
    dTdD1(i)=mean(normcdf(param(10)+param(22)+(param(14)+param(24))*theta1(i)+(param(12)+param(26))*G)...
        -normcdf(param(10)+param(14)*theta1(i)+param(12)*G));
    dTdD2(i)=mean(normcdf(param(11)+param(23)+(param(15)+param(25))*theta1(i)+(param(13)+param(27))*G)...
        -normcdf(param(11)+param(15)*theta1(i)+param(13)*G));
    JTd1 = NumJacobian(@(par)...
        mean(normcdf(par(10)+par(22)+(par(14)+par(24))*theta1(i)+(par(12)+par(26))*G)...
        -normcdf(par(10)+par(14)*theta1(i)+par(12)*G)),param);
    JTd2 = NumJacobian(@(par)...
        mean(normcdf(par(11)+par(23)+(par(15)+par(25))*theta1(i)+(par(13)+par(27))*G)...
        -normcdf(par(11)+par(15)*theta1(i)+par(13)*G)),param);
    Vtd1(i) = sqrt(JTd1*SIG*JTd1');
    Vtd2(i) = sqrt(JTd2*SIG*JTd2');
        
end    
   
[dTdD1;Vtd1;dTdD2;Vtd2]
    

    
% 
% for j=1:4;
% for i=1:100;
%     dYdTEw(j,i)=dYdT(j*2-1,1,i);
%     dUw(j,i) = dYdT(j*2-1,1,i)+1.64*dYdT(j*2,1,i);
%     dDw(j,i) = dYdT(j*2-1,1,i)-1.64*dYdT(j*2,1,i);
%     
%     dYdTEb(j,i)=dYdT(j*2-1,2,i);
%     dUb(j,i) = dYdT(j*2-1,2,i)+1.64*dYdT(j*2,2,i);
%     dDb(j,i) = dYdT(j*2-1,2,i)-1.64*dYdT(j*2,2,i);
%     
% end
% end
% 
% 
% 
% figure
% plot(theta,dYdTEb(3,:),'k',theta,dUb(3,:),'k--',theta,dDb(3,:),'k--')
% xlabel('\theta')
% ylabel('Partial Effects')
% title('PE of D_E on T_E, Black Students')
% saveas(gcf,'pedete.jpg')
% 
% figure
% plot(theta,dYdTEb(4,:),'k',theta,dUb(4,:),'k--',theta,dDb(4,:),'k--')
% xlabel('\theta')
% ylabel('Partial Effects')
% title('PE of D_M on T_M, Black Students')
% saveas(gcf,'pedmtm.jpg')
% 
% figure
% plot(theta,dYdTEw(3,:),'k',theta,dUw(3,:),'k--',theta,dDw(3,:),'k--')
% xlabel('\theta')
% ylabel('Partial Effects')
% title('PE of D_E on T_E, White Students')
% saveas(gcf,'pedetew.jpg')
% 
% figure
% plot(theta,dYdTEw(4,:),'k',theta,dUw(4,:),'k--',theta,dDw(4,:),'k--')
% xlabel('\theta')
% ylabel('Partial Effects')
% title('PE of D_M on T_M, White Students')
% saveas(gcf,'pedmtmw.jpg')
% 
% 
% 
 